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Abstract We construct several variational integrators — integrators based on a discrete vari- 
ational principle — for systems with Lagrangians of the form L = La + eLb, with e <C 1, 
where describes an integrable system. These integrators exploit that e <?C 1 to increase 
their accuracy by constructing discrete Lagrangians based on the assumption that the in- 
tegrator trajectory is close to that of the integrable system. Several of the integrators we 
present are equivalent to well-known symplectic integrators for the equivalent perturbed 
Hamiltonian systems, but their construction and error analysis is significantly simpler in the 
variational framework. One novel method we present, involving a weighted time-averaging 
of the perturbing terms, removes all errors from the integration at 6 (e). This last method is 
implicit, and involves evaluating a potentially expensive time-integral, but for some systems 
and some error tolerances it can significantly outperform traditional simulation methods. 
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1 Introduction 

Symplectic integrators have been used since their introduction in IWisdom and Holmanl ( fl99lh 
for simulations of gravitational systems which are dominated by a central body. These inte- 
grators split the Hamiltonian for a system into two parts: 

H = H {A) +£H^ B \ (1) 

where represents the influence of the dominant central body on the bodies that orbit it, 
and eH( b \ e<1, represents the mutual interactions of the bodies in orbit around it. (In our 
solar system, e ~ 10~ 3 ; for stars in orbit around a central galactic black hole, e ~ 10~ 6 .) 
These integrators are a composition of evolutions under t he separate pieces of the Ham ilto- 
nian, which are individually integrable. The integrators in lwisdom and Holmanl have 
a trajectory error which scales as & (eh ) over a single step of size /z. lMcLachlanl ( ll995h . 
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IChambers and MurisorJ §000), and Laskar and Robutel (2001) present improvements to the 



basic leapfrog scheme in IWisdom and Holmanl ([1991) which involve more composition 



steps to eliminate error terms at <?(£); these integrators have errors after a single step of 
size h which scale as 6 ( e 2 ^) + & (e// +1 ) and are known as pseudo-high-order integra- 
tors. IWisdom et al.l i ll 9961) introduced correctors which can completely eliminate errors at 
0(e) from the integration. 

In this work, we present integrators derived from a Lagrangian of the form 

L = L^+eLW (2) 

based on a discrete variational principle which incorporates the dominant motion of the 
system. These variational integrators subsume tradition al symplectic integ r ators. Section[2] 
provdies a brief introduction to variational integrators; iMarsden and West! 1 120011) provides 
a thor ough mathematica l grounding for these integrators and discusses in detail their prop- 
erties. iLew et al.l < [2004h provides another introduction to the topic, and di scusses the use 
of variational integrators in a space-time (PDE) context. iLee et al I d2007t) demonstrate a 
geometrically exact method for simulating full-body dynamics in orbital mechanics with a 
variational integrator. 

In this paper, we der i ve the pseu do-high-order integrators pr esented in lMcLachlanI dl995l) . 
IChambers and MurisorJ d2000h . and Laskar and Robutell f200lh using the variational frame- 
work. The occurrence of the Gauss-Lobatto quadrature coefficients in the composition for- 
mulas for these integrators is a natural consequence in this framework of using Gauss- 
Lobatto quadrature to approximate the contribution of to the action integral. We exploit 
the flexibility of the variational framework to derive a novel implicit integrator which elimi- 
nates all errors from the trajectory at &(b) through averaging of the perturbing Lagrangian, 
eLS s >, along the trajectory of lS A \ We present numerical evidence in Section|4]that this lat- 
ter integrator can be more efficient than standard symplectic integrators for some systems of 
physical interest in spite of its higher cost per step. We argue in Section [3~2l fhat this method 
should be more stable at large stepsize for eccentric systems than the standard symplectic 
integrators, even with symplectic correctors. 

In this work we use the notation of Sussman et al. I l l200ll) . The function which is the 
derivative of a function / is denoted Df; if / takes a vector argument x, then Df(x) is a 
co-vector. Similarly, we denote the function which is the partial derivative with respect to 
the ;th argument (counting from zero) of a function, g, of multiple arguments by dig. Some 
examples relating our notation to traditional notation: 

D/(x) = ^)or& (3) 



dx dx 

2 f , n df(x,b,z) df(x,b,z) 
<?t/( W ) = — gj— ^ or — — 



(4) 

b=y 



2 Variational Integrators 

We can construct a variational integrator for a mechanical system with Lagrangian L(g, v) — 
here assumed to be time-indepen dent for simplicity — by considering the action for the sys- 
tem over a small interval of time ( iFarr an d Bertschinger 2007;; [Lew et alj|2~004) : 



rh 

S(h,q ,qi) = dtL(q(t),Dq(t)), (5) 
Jo 
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where q(t) is the physical trajectory of the system with q(0) = qo and q(h) = q\ . We choose 
a function, called the discrete Lagrangian, which approximates the action integral: 

Lo(h,qo,qi) ~ S(h,q ,qi) . (6) 
Then the discrete Elder-Lagrange equations 

diL D (h,q ,qi) = -po (7) 
d 2 L D (h,q ,q\) = p x (8) 

define an integrator, (qo,po) i— > (q\,p\). Depending on the structure of Lp, equation l|7]l may 
be implicit for q\. 

For example, consider the Harmonic oscillator, which has Lagrangian 

L=\{v 2 -q 2 ) (9) 

in suitable units. Adopt the ansatz that q(t) = qo + (qi — qo)(t — to)/h; we can construct a 
discrete Lagrangian using the midpoint rule: 

L D (h,q ,qi)=hL\ — - — , — - — I. (10) 

The integrator equations Q and <[Sj can be solved explicitly for q\ and p\. The result is 

4-h 2 Ah 



Ah A-h 2 

Pi = -4Th- 2qo + ATi? po - (12) 

Comparing with the actual solution 

q(t) = qocos(t)+posin(t) (13) 
p(t) = -qosm{t)+p Q cos(t), (14) 

we see that the integrator follows the exact trajectory but with phase error — that is, q\ = 
q(h + 8t) and p\ = p(h + 8t) for some phase error 8t. 

One way to understand variational integrators is to recall that the action is an F\ -type 
gene rating function for the canonical transformation that implements time-evolution (see, 
e.g., ISussman et all d200ll . pp. 415—416)). The action defines the Fj-type map, (q,p) i— > 
(2,F),via 

diS{h,q,Q) = -p (15) 
d 2 S(h,q,Q) = P; (16) 

these are just equations ^ and ®, with the discrete Lagrangian — an approximate action — 
replaced by the true action. 

Alternately, consider our approximation to the action over a longer interval: 



Stot(?0,2l,?2,---) = S(h,q ,q l )+S(h,q l ,q 2 ) + ... 

w L D (h,qo,qi) +L D (h,qi,q 2 ) + . . ■ ■ (17) 
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Then equations ((TJl and $Q express the stationarity of our approximation to the action with 
respect to the intermediate positions q\,q2,.... 

<5lStot(<70,<7l,<72,---) « d 2 L D (h : q ,q l )+diL D (h,qi,q 2 ) = pi -p\ = 0. (18) 

Because the mapping in equations (fTJl and ^ is the extremization of an approximation to 
the action, it shares many of the desirable properties of the exact trajectory which extremizes 
the true action. For example, if Lq has symmetries, a discrete Nother's theorem implies 
that the corresponding momenta are conserved. Also, the mapping is symplectic: denote 
the mapping by F. Then the symplectic form on phase space dq A dp is invariant under 
pushforward by F : 

F* (dq A dp) = dq A dp. (19) 

Finally, it is possible to show via backward error analysis dLew et al.l2004) that the mapping 
implements the exact evolution for some Lagrangian L which is close to L. Therefore, the 
evolution under Lo remains on a constant-energy surface in phase space for L. Because L is 
close to L, the evolution under Lo always remains close to the constant-energy surface of L, 
so the Jong-term energy error of the mapping under Lo is bounded. 

The order of error in the mappin g of equ ations (01 and ((Hi is the same as the order of 
error in the action approximation Lo l lMarsden and WeslfeOOlT Theorem 2.3.1). That is, if 

L D (h,q(0), q{h)) =S(h,q(0), q{h)) + ff(h n+l ) , (20) 

where q(t) is a stationary-action trajectory, then the mapping defined by equations (|7]l and 
<[8j> approximates the physical trajectory to ff(h n+i ^, and therefore defines an n-th order 
integrator for L. We shall exploit this result in the error analysis of the integrators presented 
in this paper. 



2. 1 Multi-Point Variational Integrators 

It is often advantageous to allow the discrete trajectory to pass through intermediate points 
between qg and 171 . For example, we may imagine that the discrete trajectory is a polynomial 
in time which interpolates between the positions qo, q' , q'$, q^\ q\. With intermediate 



positions, the discrete Euler-Lagrange equations become 

diL D (h,q ,q' ,qQ,...,q^ ) ,qi\ = -p (21) 

di+\L D (h, q , q' , q'g, ...,q^\qi\ = 0, ?'=l,2,...,n (22) 

d„+2L D (h,q ,q' Q , q'o, . . . ,q^' ,qA = p\. (23) 



The intermediate equations express that the action is stationary with respect to the interme- 
diate positions, while the other equations give the time-evolution canonical transformation. 

We can always (in principle) re-express any multi-point discrete Lagrangian as an equiv- 
alent two-point discrete Lagrangian by solving the intermediate equations for the interme- 
diate positions in terms of the start and end positions: 

L D (h,q ,qi) =L D (h,qo 7 q' {h,q ,qi),...,q^' (h,q ,qi) ,qi) , (24) 
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where the functions q' Q (h,qo,qi), (h,qo,qi) solve 

d i+1 L D (h,q Q ,q' Q , q { Q "\q^) = 0, i=l, (25) 

Applying equations <[^0> an d © to Lo gives the same trajectory as the discrete Euler-Lagrange 
equations forZ-o. 

One way to form multi-point discrete Lagrangians is to chain together two separate 
discrete Lagrangians: 

L D (h,q ,q' ,qi) =1$ (h,q Q ,q' Q ) +L { ^ (h,q' ,q\) . (26) 

A quick calculation shows that the evolution under Lp is a composition of evolution under 
L$ with that under [Jjy . 

A common way to construct variational integrators of various orders is to assume that 
the trajectory of the system is a polynomial in time which passes through some discretization 
points: q = q (t;qo,q' , . . . ,q^\q\ \ . One then forms a discrete Lagrangian using a quadra- 
ture rule on the action integral evaluated on the discrete trajectory: 

L D (h,qo,q' Q ,...,q { "\qi^ 

= h £ W t L (a(t i \q ,q' Q ,...,q^\qA, d q (fa qo, q'o> ■ ■ ■ , ^ , 9 1 ) ) » ( 27 ) 

i 

where {w,,?,} are the weights and times of a quadrature rule. The order of the resulting map 
is determined by the orders of the polynomial interpolation and quadrature rule. 



3 Almost-Integrable Systems 

An almost-integrable system has a Lagrangian of the form 

L = L^ + eL^ B \ (28) 

where e <C 1, and the trajectories of iS^ are calculable analytically (or at least efficiently). 
In contrast to the general Lagrangian, where little can be said about trajectories, we expect 
that the trajectories of L are going to be, in some sense, close to those of L^ A \ We should 
take advantage of this fact when designing variational maps to approximate the trajectories 
of L instead of blindly assuming polynomial-in-time motion as discussed at the end of the 
last section. 

The integrators we are about to describe all use as a component a particular discrete 
Lagrangian: 

4 A) (h,q A (0),q A (h)) = J dtlW (q A (t),Dq A (t)), (29) 

where q A (t) is the trajectory corresponding to the Lagrangian Lp is the exact action 
for the A-subsystem on its trajectories; applying equations l|7]l and <[8j) to ' gives the q A 
trajectory: 

d x L { ^ (h,q A (0),q A (h)) = -p A (0) (30) 
d 2 L^ (h,q A (0),q A (h)) = PA {h). (31) 

(A) 

In general, it is not necessary to compute L D ; all that is necessary is to be able to efficiently 
compute the mapping defined by equations d30t and ( 13 1 b . 
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3.1 Composition Maps Using Quadrature Rules 

In this subsection we will show how some well-known symplectic maps are equivalent to 
variational integrators and we will analyze their errors in the variational framework. Both 
the derivation and error analysis are considerably simpler in the variational framework. We 
will assume that 

L^( q ,v) = -V^( q y, (32) 

this common case is required for the integrators we discuss to be compositiona l. 

W e will see that the well-known sym plectic integrators in lMcLachladdl995h . lChambers and MurisorJ 
J2000h . and lLaskar and Robutell fcOO lh result from using Gauss-Lobatto quadrature to ap- 
proximate 

h dM B \q{t)) (33) 



in the d iscrete Lagrangian, assum ing qa traje ctories between the quadrat ure points. In M cLachlan 
dl995h . lChambers and MurisorJ J2000h . and lLaskar and Robutell J200lh the coefficients for 
these mapping integrators come from the solution to algebraic equations involving iterated 
commutators of the Hamiltonians and , and turn out to be Gauss-Lobatto quadra- 
ture weights and times. In the framework of this paper, the Gauss-Lobatto quadrature coef- 
ficients arise naturally from the attempt to approximate the time-integral of V^ B > to a given 
order. 

3.1.1 Kick-Drift-Kick Leapfrog 
Suppose we choose 

L D (h,q ,qi) = L$ (h,q ,qi) + eL^' (h,q ,qi) 



= 4 A) (h,qo, qi )-e^ [yW (q ) + V^ (q 1 )} . (34) 
Here we have simply used the trapezoidal rule (a two-point Gauss-Lobatto quadrature) 

£dxf(x)*i^\f(a)+f{b)] (35) 

to approximate the contribution to the action from Lb- Applying equations ^ and ([Sj) to this 
discrete Lagrangian, we have 

-po = di4 A) (Mo,4i)-4w (B) (<7o) (36) 

pi =d 2 L ( g ) (h,q ,q 1 )-e^DvW(q 1 ). (37) 



These can be re- written in a suggestive way 

h 

: 2 



p -e%DVW ( qo ) ) = diL^ (h,q ,qi) (38) 



Pi =d 2 L<g ) (h,q ,qi)-e^DVW(q 1 ). (39) 



The solution to the first equation is the q\ that results from evolving yqo,po — eh/2DV^ (qo) 
forward by L~ A \ The final momentum is then p& — EjDV^ (qi). In other words, we kick 
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by -e^DVW (q ), drift by Z.M, and then kick by -e^D V^ («?,). The algorithm is kick- 
drift-kick leapfrog — see, e.g. I Wisdom and Holmanl dl99 lh . 

Using the variational framework to analyze the error of kick-drift-kick leapfrog, we 
consider 

A=L D (h,q(0),q(h))-S{h,q{0),q(h)). (40) 

The trajectory error of kick-drift-kick leapfrog will be of the same order as A . Expanding, 
we have 



A =A A + £A B 



dt [L A (q A (t;q(0),q(h)),Dq A (t;q(0),q(h)))-L A (q(t),Dq(t))} 



\ (vW(q(0))+vW(q(h))) -J H dtVW{q{t)) 



(41) 



For the first term, we have 



A A - 



8q A 



+ 



dt L A (q A (t;q(0),q(h)),Dq A (t;q(0),q{h))) 

h 



Sq A 



[ dtL A (q A {f,q(0),q(h)),Dq A (t;q(0),q{h))) 
Jo 



8q 2 A + 0(5q 3 A ), (42) 



where 8q A is the trajectory difference between q A and q. 8q A is ff (eh) because the La- 
grangian depends on the velocity, and the two trajectories must differ at first order in h in 
their velocities (since they feel different forces at order e). But, the trajectory q A is the so- 
lution to the Euler-Lagrange equations for L A , so the first order variation of L A vanishes on 
q A , and only the second-order term contributes. Therefore, we have 



A A = <ff(e 2 h 3 ). 
For the second term of equation d4 1 b . we have 

£A B = e[<ff(h 3 )} =ff(eh 3 ) 



(43) 



(44) 



arising from the truncation error in the quadrature rule. Putting the two terms together, we 
see that 



A = A A + eA B = (eh 3 ) + (e 2 h 3 ) , 



(45) 



with the term at ff (e) arising from the quadrature rule error, and the term at & (e 2 ) arising 

(A) 

from the error in L D ' . This will be a general feature of the integrators in this section: the 
quadrature error determines the ff (e) integrator error, while the (e 2 ) error is determined 

(A) 

by the error in the L D term. 
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3.1.2 S4B 

Consider choosing a higher-order quadrature rule for the eL^ term. For example, a three- 
point Gauss-Lobatto quadrature rule, with q& trajectories between the quadrature points: 

L D (h,q ,qf ,qi) = 

1$ Q>Wo) +1$ -4 [y {B) M+4vW (^ )+vW( ? i)] (46) 

As above, this integrator can be written as a sequence of kicks and drifts: a kick by —e^DV (qo), 
a drift by h/2 with respect to L' A ' to ^q, a kick by — e^DV (q' ), a drift by h/2 with respect 
toZ/ A ) to ^1, and a final kick by -e%DV(qi). 

As before the error introduced by theL^ -* part is second order in 8qA, or (e 2 h 3 ), while 
the quadrature rule introduces an error in theVW part of iff (eh 5 ). (Th ere is an additional 
error in the quadrature arising from the trajectory error SqA, which contributes at (e 2 /i 3 ) 
because the potential is independent of the velocity difference between q and q&.) 
Thus, the error in this method scales as 

eh^+0(e 2 h i ). (47) 

It belo ngs to a class of integrators known as "pseudo-high-or der", first discovered bvlMcLachlanl 
dl995), and introduced t o the astronomical community in Chambers and Murisonl d200of) 



and lLaskar and Robutell < f200lh . These integrators are useful because, for small enough e, 
they behave as high-order integrators, even thou gh they are formally second ord er. The name 
of t his section ( a nd the integrator) is S4B, from lchambers and Murisonl bOOOl) . 

iMcLachlanl (l995) originally derived the coefficients of the drifts and kicks for this 
integrator by attempting to eliminate commutator terms in the Lie series for Hamiltonian 
evolution; he noted that the coefficients which eliminate the desired first-order-in-e terms 
are identical to the Gauss-Lobatto quadrature coefficients. In this work, we can understand 
this as a consequence of attempting high-order quadrature of the contribution of to the 
action. 

3.1.3 S6B 

Consider now the sixth-order Gauss-Lobatto quadrature for the iffl terms, interspersed with 
qA evolution: 

L D (h,qo,q l Q,q l Q,q\) = 

D 1 10 ^ >4Q\ +L D y-7^^0>%)+ L D I ]7j >?0)?1 I 

- £ A [v^(q )+5V^(q' )+5V^(q'i)+vW(q 1 )] (48) 

This i s exactly the sequence of drifts and kicks for the S6B integrator ^Chambers and Murisonl 
2000). Once again, the error is a combination of errors from the quadrature at 0(s), and er- 
rors from 1$ at 6 (e 2 ) : 

A = 0(eh n ) + 0(e 2 h 3 ) (49) 



9 



3.2 An & (e 2 ) Method 

We can eliminate all errors at ff (e) by using for Lfi' the exact integral of along the q A 
trajectory. Define 

L D (h, qo,qi)= L%' (h,q ,qi) + eL (B) (h, qo,qi) 

rh 

= L<g\h,qa, qi )-e dtV^(q A (f,q ,qi)). (50) 
Jo 

Applying equation ^ to Lq, and moving the L^f' term to the left-hand-side, we obtain 

- (po + ediL^ (h,qo,qiU = L$ (h,q ,qi) , (51) 
which we must solve for q\ . The momentum kick, 

rh 

ed l L i p ) (h,q ,qi) = -e dtDV^ (q A {t;q ,qi))diq A [t;q ,qi) , (52) 
Jo 

is the time-averaged force along an LS A > trajectory weighted by d\q A (t;qo,qi) — in general, 
this weight favors the initial periods of the trajectory, since q A (h;qo,qi) = q\ independent 
of qo. Once the point q\ is determined, the new momentum is 

Pi = d 2 L D (h,q ,qi) = + ed 2 L^ (h,q Q ,qi) . (53) 

Here, the kick 

rh 

ed 2 L { * ) (h,q ,qi) = -£ dtDV (B) (q A (t;q ,qi))d 2 q A (t;q ,q 1 ) , (54) 
Jo 

is the time-averaged force along an lS A > trajectory weighted by d 2 q A (t; qo,q\ ) — which tends 
to favor later points in the trajectory, since q A (0;qQ,qi) = qo independent of q\. Because the 
momentum kicks are related to the time-averaged force along the integrator trajectory, the 
integrator has a flavor of averaging. Timesteps with this integrator can be as large as the time 
it takes the trajectory to deviate on average from q A , in contrast to the integrators from the 
previous subsection. In those integrators timesteps must be small enough that both q A ade- 
quately approximates the trajectory and that the sequence of kicks adequately approximates 
the averaged force. Evaluating the averaged force in the way that the € (e 2 ) variational 
method does removes this second restriction. This could be important in the simulation of 
highly eccentric systems. 

3.2.1 Error Analysis 

The error from the integral is 

8q A + ff(8q 2 A ) =&(e 2 h 3 ). (55) 

Because depends only on q and not on q, 8q A scales as ff (e^ 2 ) (the true trajectory q 
and q A must differ at order eh 2 because they feel different forces of size e). 



8q A 



dtVW (q A (f,q ,qi)) 
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Combining the error in equation d55l l with Aa, we obtain 

A = ff{e 2 h i ), (56) 

resulting in a method which is formally second-order, but has no errors at is.). 
The method is implicit; one must solve equation Q, 

- po = diL D (h,qo,qi) , (57) 

for q\ . This can be accomplished through Newton iteration, or through the iterative method 
which follows. (In practice, for small e < 10~ 3 and systems of modest dimension, we find 
that the iterative method is more efficient than Newton iteration.) Define the sequence |^ j 
by 

- (po-ediL^ (li,q ,q [ ;- l) )) = dil$ (h,q ,qf) . (58) 
If<7[' is known, then q^ is just the evolution of the state (qo,po — ed\l^ (h,qo,q( 

by La- For small e, this sequence (q^f J converges linearly to the desired q\ . 
The efficiency of this method will depend on how tractable it is to evaluate 

[ h dtV {B '>(qA{t;qo,qi)) (59) 

Jo 

as a function of the endpoints qo and q\. In the next section, we will examine the perfor- 
mance of the methods introduced in this section on some example problems. 

4 Example Calculations 

In this section, we apply the integrators from the previous section to some example prob- 
lems. 



4.1 The Perturbed SHO 
Consider the Lagrangian 

L(q,v) = ^{v 2 -q 2 )-^q\ (60) 

This represents a simple harmonic oscillator (with natural period 2k) with an additional 
force F(q) = —eq 2 . For £ <gC 1, the system is amenable to solution using the methods from 
the previous section. In particular, because the perturbation term is a polynomial in q, we 
can easily compute the ff (e 2 ) discrete Lagrangian: 



rh 

L D (h,q ,qi) = / dt L(q A (f,qo,qi) ,Dq A (t;q Q ,qi)) 
Jo 



/o 

~ (<7o + ?1 ) COt(ft) - q qi csc(h) 



2 

- e [(qo + qi ) (2q 2 + q qi + 2q 2 ) + (gg + q\) COS (ft)] sec 2 (j) tan (^) (61) 

Figure [T] displays the maximum energy error over a simulation of the oscillator with a 
total time T = 1000 as a function of the timestep h for the various methods in Section[3] We 
can see in Figure Q] that the ff (e 2 ) variational method significantly outperforms the other 
methods at large timesteps. 
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Fig. 1 Maximum relative energy error in a simulation of the perturbed SHO over a total time T = 1000 as 
a function of timestep, h by the methods of SectionfJ] The data are for KDK (plus), S4B (circle), S6B (star), 
and (e 2 ) variational (dot). We set e = 10~ 5 . Though it is an implicit method, and takes more work per step, 
the ff (e 2 ) variational integrator has such dramatically better error behavior at large stepsize that it is more 
efficient for integrating this system. 



4.2 Jupiter, Saturn, and the Sun 

This subsection reports on simulations of the Jupiter-Saturn-Sun system with realistic initial 
conditions. The Lagrangian for this system is 

\ , 2 2 2 \ , Gm e mj Gm e m s Gmjm s 

{inpjVQ+mjVj + m s v s ) H 1 1 . (62) 

2 r QJ r QS r JS 

The well-known Jacobi transformation (see, e.g. Wis dom and Holmanl ( 1 199 lb ) can trans- 
form this Lagrangian into a sum of center-of-mass motion, two Kepler Lagrangians, and 
perturbing terms with magnitude e ~ mj/rriQ ~ 10~ 3 . 

In this paper, we evaluate the time-average of the perturbing terms — which are essen- 
tially the disturbing function for the three-body problem — for the G (e 2 ) variational inte- 
grator on the qA (Kepler) trajectory using numerical quadrature. Numerical quadrature is 
adaptive; each quadrature point corresponds to a (weighted) kick at that time, and quadra- 
ture points are allocated non-uniformly on the interval [0, h] to best approximate the integral. 
With this technique, we can use the (e 2 ) variational method on highly elliptical orbits with 
a large timestep without loss of accuracy, since the quadrature routine will allocate points 
densely near pericenter. The other integration methods, which allocate quadrature points at 
fixed fractions of the stepsize, must be run with a small enough stepsize to resolve rapidly 
changing forces near pericenter passage throughout the entire orbit. 

Figure [2] displays the energy error in a simulation of the Jupiter-Saturn-Sun system 
for approximately 20 Jupiter orbits (which corresponds to approximately 240 years) ver- 
sus timestep. For a maximum tolerable error of e 2 ~ 10~ 6 , the ff (e 2 ) variational integrator 
can take stepsizes of order 10 orbits, while the other methods do not perform well until 
there are several kicks per orbit. However, for high-accuracy integrations the errors in the 
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Fig. 2 Maximum relative energy error versus stepsize in a simulation of Jupiter, Saturn and the Sun's mutual 
gravitational interaction over a period of 250 years (Jupiter's period is about 12 years). The curves are for the 
following algorithms: KDK (plus), S4B (circles), S6B (stars), and ff (e 2 ) variational (dots). Because in this 
system e ~ 2£. K 10 -3 is relatively large, the C (e 2 ) variational integrator is only advantageous at energy 

errors of order e 2 ~ 10~ 6 . However, for such large energy en'ors, the (e 2 ) variational integrator can take 
stepsizes which are on the order of 10 orbital periods, significantly larger than traditional algorithms; these 
large stepsizes more than offset the increased computational cost of the method. 



methods (excepting KDK) are comparable, and the extra cost of the averaging in the ff (e 2 ) 
variational integrator when compared to the other methods makes it sub-optimal. 



4.3 Small-mass Jupiter, Saturn, and the Sun 



This subsection reports on a simulation with the same initial conditions as Section [4721 but 
with the masses of Jupiter and Saturn reduced by a factor of 10~ 3 . This brings e ~ 10~ 6 , 
roughly in line with the size of the perturbing interaction one might find in a cluster of stars 
around a super-massive black hole in the center of a galaxy. 

Figure [3] presents the relative energy error versus timestep for a simulation of this 
smaller-e system over approximately 100 Jupiter orbits. In this circumstance, the & (e 2 ) 
variational integrator significantly outperforms the other integrators, even for the (relatively 
severe) error budget of 10~ 12 . It can take steps which are approximately 10 3 longer than 
those of the other integrators at the same energy error budget, more than compensating for 
the expensive time-averaging and implicit nature of the algorithm. 

In Figure |4] we plot the trajectory error of the various methods at the end of the the 
simulation period (100 Jupiter orbits). The ff (e 2 ) variational method outperforms the other 
methods by approximately a factor of 10 3 in stepsize at a relative trajectory error of 10~ . 
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Fig. 3 Maximum relative energy error in a simulation analogous to the one in Section l4~2l except with Jupiter 
and Saturn's masses reduced by a factor of 10~ 3 (note that Jupiter still has a 12-year period). The curves are 
for the following algorithms: KDK (plus), S4B (circles), S6B (stars), and (e 2 ) variational (dots). In this 
system e ~ 10~ 6 , and we see that the (e 2 ) variational integrator can take stepsizes which are ~ 10 3 larg er 
than other algorithms for an error budget of 10 -12 . 
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Fig. 4 Relative phase-space (trajectory) error at the end of a simulation analogous to the one in Section 
14.21 except with Jupiter and Saturn's masses reduced by a factor of 10~ 3 (note that Jupiter still has a 12- 
year period). The curves are for the following algorithms: KDK (plus), S4B (circles), S6B (stars), and (e 2 ) 
variational (dots). In this system e ~ 1(P 6 , and we see that the (e 2 ) variational integrator can take stepsizes 
which are ~ 10 3 larger than other algorithms for a phase-space error budget of 10 -10 . 
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5 Conclusion 

The variational framework subsumes standard symplectic methods. In this work, we have 
presen ted t he pseudo-high-order integ rators of lMcLachlanl l l 19951) . IChambers and Murisonl 
(2000), and lLaskar and Robutel d2001l) from the variational viewpoint for systems with La- 
grangian L = Z/ A ' + eZ/ B ' . In addition, we have used the variational framework to derive a 
novel implicit integrator which uses the average perturbing Lagrangian over trajectories of 
the dominant Lagrangian to remove all errors from the integration at ff(s). We have pre- 
sented numerical evidence that, for small e, this latter integrator is more efficient than stan- 
dard pseudo-high-order symplectic integrators for perturbed systems. It would be interesting 
to investigate the performance of this latter integrator with various analytical approximations 
to the average of the perturbing Lagrangian. 
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